Gradient Descent Methods

$$\newcommand{\NN}{\mathbb{N}} \newcommand{\CC}{\mathbb{C}} \newcommand{\GG}{\mathbb{G}} \newcommand{\LL}{\mathbb{L}} \newcommand{\PP}{\mathbb{P}} \newcommand{\QQ}{\mathbb{Q}} \newcommand{\RR}{\mathbb{R}} \newcommand{\VV}{\mathbb{V}} \newcommand{\ZZ}{\mathbb{Z}} \newcommand{\FF}{\mathbb{F}} \newcommand{\KK}{\mathbb{K}} \newcommand{\UU}{\mathbb{U}} \newcommand{\EE}{\mathbb{E}} \newcommand{\Aa}{\mathcal{A}} \newcommand{\Bb}{\mathcal{B}} \newcommand{\Cc}{\mathcal{C}} \newcommand{\Dd}{\mathcal{D}} \newcommand{\Ee}{\mathcal{E}} \newcommand{\Ff}{\mathcal{F}} \newcommand{\Gg}{\mathcal{G}} \newcommand{\Hh}{\mathcal{H}} \newcommand{\Ii}{\mathcal{I}} \newcommand{\Jj}{\mathcal{J}} \newcommand{\Kk}{\mathcal{K}} \newcommand{\Ll}{\mathcal{L}} \newcommand{\Mm}{\mathcal{M}} \newcommand{\Nn}{\mathcal{N}} \newcommand{\Oo}{\mathcal{O}} \newcommand{\Pp}{\mathcal{P}} \newcommand{\Qq}{\mathcal{Q}} \newcommand{\Rr}{\mathcal{R}} \newcommand{\Ss}{\mathcal{S}} \newcommand{\Tt}{\mathcal{T}} \newcommand{\Uu}{\mathcal{U}} \newcommand{\Vv}{\mathcal{V}} \newcommand{\Ww}{\mathcal{W}} \newcommand{\Xx}{\mathcal{X}} \newcommand{\Yy}{\mathcal{Y}} \newcommand{\Zz}{\mathcal{Z}} \newcommand{\al}{\alpha} \newcommand{\la}{\lambda} \newcommand{\ga}{\gamma} \newcommand{\Ga}{\Gamma} \newcommand{\La}{\Lambda} \newcommand{\Si}{\Sigma} \newcommand{\si}{\sigma} \newcommand{\be}{\beta} \newcommand{\de}{\delta} \newcommand{\De}{\Delta} \renewcommand{\phi}{\varphi} \renewcommand{\th}{\theta} \newcommand{\om}{\omega} \newcommand{\Om}{\Omega} \renewcommand{\epsilon}{\varepsilon} \newcommand{\Calpha}{\mathrm{C}^\al} \newcommand{\Cbeta}{\mathrm{C}^\be} \newcommand{\Cal}{\text{C}^\al} \newcommand{\Cdeux}{\text{C}^{2}} \newcommand{\Cun}{\text{C}^{1}} \newcommand{\Calt}[1]{\text{C}^{#1}} \newcommand{\lun}{\ell^1} \newcommand{\ldeux}{\ell^2} \newcommand{\linf}{\ell^\infty} \newcommand{\ldeuxj}{{\ldeux_j}} \newcommand{\Lun}{\text{\upshape L}^1} \newcommand{\Ldeux}{\text{\upshape L}^2} \newcommand{\Lp}{\text{\upshape L}^p} \newcommand{\Lq}{\text{\upshape L}^q} \newcommand{\Linf}{\text{\upshape L}^\infty} \newcommand{\lzero}{\ell^0} \newcommand{\lp}{\ell^p} \renewcommand{\d}{\ins{d}} \newcommand{\Grad}{\text{Grad}} \newcommand{\grad}{\text{grad}} \renewcommand{\div}{\text{div}} \newcommand{\diag}{\text{diag}} \newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} } \newcommand{\pdd}[2]{ \frac{ \partial^2 #1}{\partial #2^2} } \newcommand{\dotp}[2]{\langle #1,\,#2\rangle} \newcommand{\norm}[1]{\| #1 \|} \newcommand{\normi}[1]{\norm{#1}_{\infty}} \newcommand{\normu}[1]{\norm{#1}_{1}} \newcommand{\normz}[1]{\norm{#1}_{0}} \newcommand{\abs}[1]{\vert #1 \vert} \newcommand{\argmin}{\text{argmin}} \newcommand{\argmax}{\text{argmax}} \newcommand{\uargmin}[1]{\underset{#1}{\argmin}\;} \newcommand{\uargmax}[1]{\underset{#1}{\argmax}\;} \newcommand{\umin}[1]{\underset{#1}{\min}\;} \newcommand{\umax}[1]{\underset{#1}{\max}\;} \newcommand{\pa}[1]{\left( #1 \right)} \newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. } \newcommand{\enscond}[2]{ \left\{ #1 \;:\; #2 \right\} } \newcommand{\qandq}{ \quad \text{and} \quad } \newcommand{\qqandqq}{ \qquad \text{and} \qquad } \newcommand{\qifq}{ \quad \text{if} \quad } \newcommand{\qqifqq}{ \qquad \text{if} \qquad } \newcommand{\qwhereq}{ \quad \text{where} \quad } \newcommand{\qqwhereqq}{ \qquad \text{where} \qquad } \newcommand{\qwithq}{ \quad \text{with} \quad } \newcommand{\qqwithqq}{ \qquad \text{with} \qquad } \newcommand{\qforq}{ \quad \text{for} \quad } \newcommand{\qqforqq}{ \qquad \text{for} \qquad } \newcommand{\qqsinceqq}{ \qquad \text{since} \qquad } \newcommand{\qsinceq}{ \quad \text{since} \quad } \newcommand{\qarrq}{\quad\Longrightarrow\quad} \newcommand{\qqarrqq}{\quad\Longrightarrow\quad} \newcommand{\qiffq}{\quad\Longleftrightarrow\quad} \newcommand{\qqiffqq}{\qquad\Longleftrightarrow\qquad} \newcommand{\qsubjq}{ \quad \text{subject to} \quad } \newcommand{\qqsubjqq}{ \qquad \text{subject to} \qquad }$$

This tour explores the use of gradient descent method for unconstrained and constrained optimization of a smooth function

Contents

Installing toolboxes and setting up the path.

You need to download the following files: signal toolbox and general toolbox.

You need to unzip these toolboxes in your working directory, so that you have toolbox_signal and toolbox_general in your directory.

For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.

Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.

Execute this line only if you are using Matlab.

getd = @(p)path(p,path); % scilab users must *not* execute this

Then you can add the toolboxes to the path.

getd('toolbox_signal/');
getd('toolbox_general/');

Gradient Descent for Unconstrained Problems

We consider the problem of finding a minimum of a function f$f$, hence solving

$$\umin{x \in \RR^d} f(x)$$ where f:Rd→R$f : \RR^d \rightarrow \RR$ is a smooth function.

Note that the minimum is not necessarily unique. In the general case, f$f$ might exhibit local minima, in which case the proposed algorithms is not expected to find a global minimizer of the problem. In this tour, we restrict our attention to convex function, so that the methods will converge to a global minimizer.

The simplest method is the gradient descent, that computes

$$x^{(k+1)} = x^{(k)} - \tau_k \nabla f(x^{(k)}),$$ where τk>0$\tau_k>0$ is a step size, and ∇f(x)∈Rd$\nabla f(x) \in \RR^d$ is the gradient of f$f$ at the point x$x$, and x(0)∈Rd$x^{(0)} \in \RR^d$ is any initial point.

In the convex case, if f$f$ is of class C2$C^2$, in order to ensure convergence, the step size should satisfy

$$0 < \tau_k < \frac{2}{ \sup_x \norm{Hf(x)} }$$ where Hf(x)∈Rd×d$Hf(x) \in \RR^{d \times d}$ is the Hessian of f$f$ at x$x$ and ∥⋅∥$\norm{\cdot}$ is the spectral operator norm (largest eigenvalue).

Gradient Descent in 2-D

We consider a simple problem, corresponding to the minimization of a 2-D quadratic form

$$f(x) = \frac{1}{2} \pa{ x_1^2 + \eta x_2^2, }$$ where η>0$\eta>0$ controls the anisotropy, and hence the difficulty, of the problem.

Anisotropy parameter η$\eta$.

eta = 10;

Function f$f$.

f = @(x)( x(1)^2 + eta*x(2)^2 ) /2;

Background image of the function.

t = linspace(-.7,.7,101);
[u,v] = meshgrid(t,t);
F = ( u.^2 + eta*v.^2 )/2 ;

Display the function as a 2-D image.

clf; hold on;
imagesc(t,t,F); colormap jet(256);
contour(t,t,F, 20, 'k');
axis off; axis equal;

Gradient.

Gradf = @(x)[x(1); eta*x(2)];

The step size should satisfy τk<2/η$\tau_k < 2/\eta$. We use here a constrant step size.

tau = 1.8/eta;

Exercice 1: (the solution is exo1.m) Perform the gradient descent using a fixed step size τk=τ$\tau_k=\tau$. Display the decay of the energy f(x(k))$f(x^{(k)})$ through the iteration. Save the iterates so that X(:,k) corresponds to x(k)$x^{(k)}$.

exo1;

Display the iterations.

clf; hold on;
imagesc(t,t,F); colormap jet(256);
contour(t,t,F, 20, 'k');
h = plot(X(1,:), X(2,:), 'k.-');
set(h, 'LineWidth', 2);
set(h, 'MarkerSize', 15);
axis off; axis equal;

Exercice 2: (the solution is exo2.m) Display the iteration for several different step sizes.

exo2;

Gradient Descent in Image Processing

We consider now the problem of denoising an image y∈Rd$y \in \RR^d$ where d=n×n$d = n \times n$ is the number of pixels (n$n$ being the number of rows/columns in the image).

Load a clean image.

n = 256;
name = 'lena';
x0 = rescale( load_image(name, n) );

Display it.

clf;
imageplot(x0);

Add noise to it, to simulate a noisy image.

sigma = .1;
y = x0 + randn(n)*sigma;

Display the noisy image y$y$.

clf;
imageplot(clamp(y));

Denoising is obtained by minimizing the following functional

$$\umin{x \in \RR^d} f(x) = \frac{1}{2} \norm{ y-x}^2 + \la J(x)$$ where J(x)$J(x)$ is a smoothed total variation of the image. $$J(x) = \sum_i \norm{ (G x)_i }_{\epsilon}$$ where (Gx)i∈R2$(Gx)_i \in \RR^2$ is an approximation of the gradient of x$x$ at pixel i$i$ and for u∈R2$u \in \RR^2$, we use the following smoothing of the L2$L^2$ norm in R2$\RR^2$ $$\norm{u}_\epsilon = \sqrt{ \epsilon^2 + \norm{u}^2 },$$ for a small value of ε>0$\epsilon>0$.

The gradient of the functional read

$$\nabla f(x) = x-y + \lambda \nabla J(x)$$ where the gradient of the smoothed TV norm is $$\nabla J(x)_i = G^*( u ) \qwhereq u_i = \frac{ (G x)_i }{\norm{ (G x)_i }_\epsilon}$$ where G∗$G^*$ is the adjoint operator of G$G$ which corresponds to minus a discretized divergence.

Value for λ$\lambda$.

lambda = .3;

Value for ε$\epsilon$.

epsilon = 1e-3;

TV norm.

Amplitude = @(u)sqrt(epsilon^2 + sum(u.^2,3));
J = @(x)sum(sum(Amplitude(Grad(x))));

Function to minimize.

f = @(x)1/2*norm(x-y,'fro')^2 + lambda*J(x);

Gradient of TV norm. Note that div implement −G∗$-G^*$.

Normalize = @(u)u./repmat(Amplitude(u), [1 1 2]);
GradJ = @(x)-div( Normalize(grad(x)) );

Gradient of the functional.

Gradf = @(x)x-y+lambda*GradJ(x);

The step size should satisfy

$$0 < \tau_k < \frac{2}{ 1 + 4 \lambda / \epsilon }.$$ Here we use a slightly larger step size, which still work in practice.

tau = 1.8/( 1 + lambda*8/epsilon );
tau = tau*4;

Exercice 3: (the solution is exo3.m) Implement the gradient descent. Monitor the decay of f$f$ through the iterations.

exo3;

Display the resulting denoised image.

clf;
imageplot(x);

Copyright (c) 2010 Gabriel Peyre